;; make topography file for Jablonowski & Williamson (2006) dynamical core test

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

fname_out="./orog_lsm_alb.t30.grd"

nlat=48
nlon=96

;; set coordinates
gau_info=gaus(nlat/2)
lat=tofloat(gau_info(::-1,0))
dlon=360./nlon
lon=dlon*ispan(0,nlon-1,1)

;; orography 
gh = new((/nlat,nlon/),"float")
gh!0="lat"
gh!1="lon"
gh&lat=lat
gh&lon=lon
;; land-sea mask
lsm = new((/nlat,nlon/),"float")
copy_VarCoords(gh,lsm)
;; albedo
alb = new((/nlat,nlon/),"float")
copy_VarCoords(gh,alb)

;; make entire globe ocean
lsm = 0
alb = 7.00226

;; surface height
Omega=7.29212e-5
a=6.371229e+6
g=9.80616
u0=35
pi=4*atan2(1,1)
eta_0=0.252
eta_v=(1-eta_0)*(pi/2)   ;; Eq.(1)
c = cos(eta_v)^1.5
phi=lat*(pi/180)
cosphi=cos(phi)
sinphi=sin(phi)
;;; Eq.(7)
zs = u0*c*( (-2*(sinphi^6)*(cosphi^2+1./3.)+10./63.)*u0*c + ((8./5.)*(cosphi^3*(sinphi^2+2./3.))-pi/4.)*a*Omega)/g

gh = conform(gh, zs, 0)

;;; output to file
setfileoption("bin","WriteByteOrder","BigEndian")
fbinrecwrite(fname_out, -1,  gh)
fbinrecwrite(fname_out, -1, lsm)
fbinrecwrite(fname_out, -1, alb)
